close all; clear all; path(pathdef); clc;

% add subfunctions
addpath(genpath('subfunctions/'))
addpath(genpath('COMSOL_simulation/'))
addpath(genpath('data_figure2/'))

% color library
color_list{1} = [67 84 147]/255;        % dark blue
color_list{2} = [110 153 201]/255;      % light blue
color_list{3} = [196 40 27]/255;        % dark red
color_list{4} = [218 134 121]/255;      % light red
color_list{5} = [40 127 70]/255;        % dark green
color_list{6} = [161 196 139]/255;      % light green
color_list{7} = [218 165 32]/255;       % gold
color_list{8} = [238 221 130]/255;      % light gold
color_list{9} = [128 0 128]/255;        % purple
color_list{10} = [204 153 255]/255;     % light purple
color_list{11} = [161 165 162]/255;     % grey

%% Load data

% load figure 2 data sets
fn = '150mA-2017-11-18_08-33-17/gResultData.mat';
load(fn)
gResultSave = gResultSave(1:end-1);
image_fn = gResultSave{1}.gSeqParameter.filename;
heat_pos = gResultSave{1}.gSeqParameter.HeatingPosition_ImageNVC/gResultSave{1}.gScan.vol_to_dist/1e6;
info_names = {'NVx','NVy','NVz','sens_fin','fl','fcen','linewidth','slope',...
    'SumTemp_on','SumTemp_on_err','SumTemp_off','SumTemp_off_err',...
    'SumTemp_diff','SumTemp_diff_err','ShotNoise_exp','ShotNoise_th'};%,'GreenmA','InfraredmA'};
file_type = 'gradient';
fontsize = 16;

info_list = load_temp_data(file_type,info_names,{gResultSave});
info_list(:,1:3) = info_list(:,1:3)/1e6;
info_list = info_list(sum((info_list(:,1:2)-repmat(heat_pos(1:2),[size(info_list,1),1])).^2,2)<2e-8,:);%what is the limit setting? 2e-10


%% Heat distribution

%%% extracted from experiment %%%
% first make plot which only takes distance into account
dist = sqrt(sum(([info_list(:,1)-heat_pos(1),info_list(:,2)-heat_pos(2)]).^2,2))-0e-6;
figure(11); clf; set(gcf,'position',[100,100,800,400])
subplot(1,3,1)
errorbar(dist*1e6, info_list(:,9), info_list(:,10), 'o')
axisinfo('Distance (\mum)', 'Ton (K)', [], [], 12)
subplot(1,3,2)
errorbar(dist*1e6, info_list(:,11), info_list(:,12), 'o')
axisinfo('Distance (\mum)', 'Toff (K)', [], [], 12)
subplot(1,3,3)
errorbar(dist*1e6, info_list(:,13), info_list(:,14), 'o')
axisinfo('Distance (\mum)', 'Tdiff (K)', [], [], 12)

%%% comparing to heating simulations %%%
P = struct(); P.heatz = 0e-5;
z_list = 10e-5;
% the simulations may have arbitrary factors associated with transmission
% at various points in the system and other parameters of the sort.
% Therefore, we fit to the profile instead.
x_sim = linspace(0,max(abs(dist)),100);
fit_err = zeros([2,length(z_list)]);
beta_err = zeros([2,length(z_list)]);
beta = zeros([2,length(z_list)]);
% optimize over different possible heights for now
for i = 1
    for j = 1:length(z_list)
        P.ins_bot = i;
        P.heatz = z_list(j);
        [data_sim,P] = load_sim_data(P);
        y = P.heatz;    % conversion into units of um
        ind = abs(data_sim(:,2)-y)<2e-7;
        % isotropic fitting
        % fit error is actually now the uncertainty
        [beta(i+1,j), beta_err(i+1,j), gof] = fit_scaling_factor(data_sim(ind,1),data_sim(ind,6),dist,info_list(:,13),info_list(:,14));
        fit_err(i+1,j) = gof.rsquare;
    end
end
[~,j] = min(min(fit_err));

P.ins_bot = 1;
P.heatz = z_list(j);
[data_sim,P] = load_sim_data(P);
y = P.heatz;
ind = abs(data_sim(:,2)-y)<2e-7;
[beta0, beta0_err] = fit_scaling_factor(data_sim(ind,1),data_sim(ind,6),dist,info_list(:,13),info_list(:,14));
T_sim = interp1(data_sim(ind,1),data_sim(ind,6).*beta0,x_sim);

% 1D temperature plot
figure(11)
subplot(1,3,3)
hold on
plot(x_sim*1e6,T_sim,'r')
title(['beta=',num2str(beta0,'%.3f'),'\pm',num2str(beta0_err,'%.3f'),...
    ',type=',num2str(i)],'fontsize',fontsize)
%% 2D heat map - experiment

new_image_fn = ['data_figure2/', fn(1:6), image_fn(40:end)];

dx = (max(info_list(:,1))-min(info_list(:,1)))/50;
[xq,yq] = meshgrid(1.5*min(info_list(:,1)):dx:1.5*max(info_list(:,1)),...
    1.5*min(info_list(:,2)):dx:1.5*max(info_list(:,2)));
vq = griddata(info_list(:,1),info_list(:,2),info_list(:,13),xq,yq,'natural');

figure(12); clf; set(gcf,'position',[300,300,600,800])
axt = subplot(2,1,1); axt2 = subplot(2,1,2);
colormap('jet')

set(axt,'units','points'); set(axt2,'units','points');
p1 = get(axt,'pos'); p1(1)=20; p2 = get(axt2,'pos'); p2(1)=20;
delete(axt); delete(axt2);
[~,ax1] = plot_confocal(new_image_fn,fontsize);
set(ax1,'units','points'); set(ax1,'pos',p1)
ax2 = axes; axis(ax2,'off'); set(ax2,'units','points');
colormap(ax2,jet);

[~,~] = contourfm(yq*1e6,xq*1e6,vq,'w','LevelStep',round((max(max(vq))-min(min(vq)))/20),...
    'LabelSpacing',1000,'ShowText','off', 'linewidth', 1);

colorbar;
caxis([8 16])
daspect([1,1,1])
set(gca, 'color', 'none')
title('Experimental result','fontsize',fontsize); set(ax2,'pos',p1)
linkaxes([ax1,ax2])
hold on
plot(heat_pos(1)*1e6,heat_pos(2)*1e6,'x','markersize',20, 'color', 'yellow', 'linewidth', 2)
%% 2D heat map - simulation
figure(12)
[h1,ax1] = plot_confocal(new_image_fn,fontsize);
set(ax1,'units','points'); set(ax1,'pos',p2)
ax2 = axes; axis(ax2,'off'); set(ax2,'units','points');
colormap(ax2,jet);
distq = sqrt((xq-heat_pos(1)).^2+(yq-heat_pos(2)).^2); Tq = beta0*interp1(data_sim(ind,1),data_sim(ind,6),distq);
[~,c2] = contourfm(yq*1e6,xq*1e6,Tq,'w','LevelStep',round((max(max(Tq))-min(min(Tq)))/20),...
    'LabelSpacing',1000,'ShowText','off', 'linewidth', 1);

colorbar
caxis([8 16])
daspect([1,1,1])
set(gca, 'color', 'none')
title('Simulated result','fontsize',fontsize); set(ax2,'pos',p2)
linkaxes([ax1,ax2])
hold on
plot(heat_pos(1)*1e6,heat_pos(2)*1e6,'x','markersize',20, 'color', 'yellow', 'linewidth', 2)
set(gca, 'color', 'none')

%% Optimization of choice of effective heating location
% optimize shift to get best fit
[~,~,err0]      = fit_scaling_factor(data_sim(ind,1), data_sim(ind,6), sqrt((info_list(:,1)-heat_pos(1)).^2+(info_list(:,2)-heat_pos(2)).^2)', info_list(:,13), info_list(:,14));
err0            = err0.rmse;
[shift, err]    = fit_shift(data_sim(ind,1), data_sim(ind,6), [info_list(:,1)-heat_pos(1),info_list(:,2)-heat_pos(2)]', info_list(:,13), info_list(:,14));

%% After optimizing heating location, compare fit result again
P.ins_bot = 1;
P.heatz = z_list(j);
[data_sim,P] = load_sim_data(P);
y = P.heatz;
ind = abs(data_sim(:,2)-y)<2e-7;
[beta0, beta0_err] = fit_scaling_factor(data_sim(ind,1),data_sim(ind,6),sqrt((info_list(:,1)-heat_pos(1)+shift(1)).^2+(info_list(:,2)-heat_pos(2)+shift(2)).^2)',info_list(:,13),info_list(:,14));
T_sim = interp1(data_sim(ind,1),data_sim(ind,6).*beta0,x_sim);

%%% 1D temperature plot
bPlotOption.MarkerSize  = 5;
bPlotOption.CapSize     = 0;
bPlotOption.LineWidth   = 1;
bPlotOption.FontSize    = 12;
bPlotOption.MarkerList = {'o', 's'};

fh = figure(15); clf;
set(fh, 'position', [100 100 [185 200]])

% simulation
plot(x_sim*1e6,T_sim, 'color', color_list{6}, 'linewidth', 2)
hold on

xdata       = sqrt((info_list(:,1)-heat_pos(1)+shift(1)).^2+(info_list(:,2)-heat_pos(2)+shift(2)).^2)'*1e6;
ydata       = info_list(:,13);
ydataerr    = info_list(:,14);

plotcustom(xdata, ydata, ydataerr, bPlotOption, 1, 5)
axisinfo('Distance (\mum)', '\DeltaT (K)', [0 21], [6 17], bPlotOption.FontSize)
set(gca, 'XTick', [0 10 20])
set(gca, 'YTick', [0 6 10 14])

%% Power scaling

% IR calibration
LD_reading                  = [0 25 50 100 150 200 250]; % mA
Power_mW_before_obj         = [0 0.92 3.85 9.92 15.9 21.9 27.8]; % mW
transmission                = 0.35;
mA_to_mW                    = fitplot({LD_reading(2:end), Power_mW_before_obj(2:end)}, 'linear', 'noplot');

Pow_vals  = transmission * mA_to_mW([50 100 150 200]);
Tmax_vals = [3.94 9.6 15.4  21.6];
Tmax_errs = [0.01 0.08 0.26 1.16];

fh = figure(16); clf;
set(fh, 'position', [100 100 [80 70]])

fitplot({Pow_vals, Tmax_vals, Tmax_errs}, 'linear', color_list{6})
hold on
plotcustom(Pow_vals, Tmax_vals, Tmax_errs, bPlotOption, 1, 5)
set(gca, 'XTick', [0 5 10])
set(gca, 'YTick', [0 10 20])
axisinfo('P (mW)', '\DeltaT_{max} (K)', [0 9], [0 25], 8)


%% Shot noise limit comparison

load('HeatGradient-2017-11-25_20-31-23/gResultData.mat')
gResultSave = gResultSave(1:end-1);
image_fn = gResultSave{1}.gSeqParameter.filename;
heat_pos = gResultSave{1}.gSeqParameter.HeatingPosition_ImageNVC/gResultSave{1}.gScan.vol_to_dist/1e6;
info_names = {'NVx','NVy','NVz','sens_fin','fl','fcen','linewidth','slope',...
    'SumTemp_on','SumTemp_on_err','SumTemp_off','SumTemp_off_err',...
    'SumTemp_diff','SumTemp_diff_err','ShotNoise_exp','ShotNoise_th'};%,'GreenmA','InfraredmA'};
file_type = 'gradient';
fontsize = 16;

info_list = load_temp_data(file_type,info_names,{gResultSave});
info_list(:,1:3) = info_list(:,1:3)/1e6;
info_list = info_list(sum((info_list(:,1:2)-repmat(heat_pos(1:2),[size(info_list,1),1])).^2,2)<2e-8,:);%what is the limit setting? 2e-10

% plot
figure(19); clf; set(gcf,'position',[100,100,[300,250]])
plot([1 10], [1 10], '-', 'color', color_list{4}, 'linewidth', 2)
hold on
plotcustom(info_list(:,15),info_list(:,16), info_list(:,16).*info_list(:,14)./info_list(:,13), [], 1, 1)
axisinfo('\eta_{shot} (K/Hz^{1/2})', '\eta_{exp} (K/Hz^{1/2})', [1 10], [1 10], 12)
set(gca, 'XTick', [1 4 7 10])
set(gca, 'YTick', [1 4 7 10])
